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Abstract 


An implicit characteristic-based approach for numerical solution of Maxwell's time-dependent curl 
equations influx conservative form is introduced. This method combines a characteristic based finite differ- 
ence spatial approximation with an implicit lower-upper approximate factorization (LU/AF) time integra- 
tion scheme. This approach is advantageous for three-dimensional applications because the characteristic 
differencing enables a two-factor approximate factorization that retains its unconditional stability in three 
space dimensions, and it does not require solution of tridiagonal systems. Results are given both for a Fourier 
analysis of stability, damping and dispersion properties, and for one-dimensional model problems involv- 
ing propagation and scattering for free space and dielectric materials using both uniform and nonuniform 
grids. The explicit FDTD algorithm is used as a convenient reference algorithm for comparison. The one- 
dimensional results indicate that for low frequency problems on a highly resolved uniform or nonuniform 
grid, this LU/AF algorithm can produce accurate solutions at Courant numbers significantly greater than 
one, with a corresponding improvement in efficiency for simulating a given period of time. This approach 
appears promising for development of LU/AF schemes for three dimensional applications. 

1 Introduction 

1.1 Background 

In the regime of finite-difference solutions for Maxwell's time dependent curl equations, ex- 
plicit methods such as the Finite-Difference Time-Domain (FDTD) method [1] — [4], the Transmis- 
sion Line Method (TLM), or (more recently) Finite- Volume Time-Domain (FVTD) methods have 
been standard techniques for the past 20-30 years. These methods are relatively simple and effi- 
cient, and have proven very robust and adequate for many types of problems. 

Since these explicit methods are conditionally stable, their maximum time step is limited by 
a stability constraint that depends on the local grid spacing. The conditional stability restriction 
can increase the cost of analyzing electromagnetic (EM) problems, especially when nonuniform 
grids are needed, because a small local grid structure can require a time step smaller than would 
otherwise be required to accurately resolve the physical wave propagation. On the other hand, 
implicit algorithms can provide unconditional stability which allows the time step to be selected 
for accuracy and resolution of the frequency content of the excitation sources, without a stability 
restriction. 

Implicit algorithms were first proposed for electromagnetics by Holland [5]. They have not re- 
ceived widespread use in electromagnetics, however, perhaps due to increased complexity and in 
some cases the need to solve a linear system of equations. For example, centered spatial difference 
approximations in conjunction with alternating direction implicit (ADI) schemes require the solu- 
tion of (for example) tridiagonal systems for each coordinate direction. Furthermore, although this 
type of scheme is unconditionally stable for the first-order wave equation in two dimensions, it is 
known to be unconditionally unstable in three dimensions [6]. The usual Fourier stability analysis 
does not account for nonperiodic boundary conditions, however, and it has been shown using a 
matrix stability analysis that this three dimensional centered difference algorithm is conditionally 


1 



stable when used with upwind boundary conditions [7]. 

There has been recent work to develop characteristic based explicit and implicit algorithms 
for CEM that are related to some solution algorithms from computational fluid dynamics. These 
methods have been developed for curvilinear three dimensional grids, using both finite difference 
and finite volume approximations. Shankar et al. [8], [9] developed explicit characteristic-based 
finite-volume methods, using both flux-vector and flux-difference splittings for electromagnetics. 
Shang, et. al. [10]— [19] have developed both implicit and explicit characteristic based methods. 
Shang [10]— [11] developed an implicit characteristic based ADI method that leads to simple bidi- 
agonal systems that are easily solved, rather than tridiagonal or pentadiagonal systems. This 
characteristic-based scheme is unconditionally stable in two dimensions and conditionally stable 
in three dimensions [11]. Shang and Fithen [18] also developed characteristic based finite differ- 
ence and finite volume schemes using Runge-Kutta time integration, in curvilinear coordinates. 
Gaitonde and Shang [20] and Turkel [4] have recently explored high order accurate compact dif- 
ferencing schemes that require tridiagonal solutions for spatial approximations, whether implicit 
or explicit time integration is used. Turkel [4] has suggested a two dimensional implicit version 
based on an ADI scheme. More recent work includes development of two and three-dimensional 
implicit ADI FDTD schemes [21]— [24]. However, these ADI schemes still require solution of a 
tridiagonal system of equations. 

1.2 Present Work 

A characteristic-based approach to spatial differencing has advantages in constructing implicit 
algorithms, as well as in utilizing the natural method of incorporating boundary conditions for 
hyperbolic systems. The ability of characteristic based methods to provide accurate nonreflecting 
solutions at outer boundaries is well known and is discussed for electromagnetics by Shang [10]- 
[11] and others. 

The present method combines a characteristic based approach to spatial differencing with an 
implicit lower- upper approximate factorization (LU/AF) time integration scheme that can be de- 
veloped as a two-factor scheme that retains its unconditional stability in three space dimensions, 
unlike three-factor ADI schemes. The LU/AF scheme also avoids the need to solve tridiagonal 
implicit systems. A two-factor approximate factorization for three dimensions was first suggested 
by Jameson and Turkel [25]. This type of factorization (two implicit passes in three dimensions) 
has been combined with characteristic based differencing by Belk and Whitfield [26]. This general 
approach has been further developed for CFD and is discussed, for example, in [27]. 

In the present report, an implicit characteristic based algorithm is developed for the time de- 
pendent Maxwell equations in one space dimension. This (one-dimensional) algorithm is investi- 
gated along two lines: A Fourier analysis is used to study both stability and dispersion/damping 
properties of the algorithm. Although the Fourier analysis provides useful information and is 
widely used to study algorithms, it is limited to periodic solutions on uniform grids. Complex 
practical applications often involve nonperiodic solutions in finite regions with complex geome- 
try and nonuniform grids. Grid nonuniformity is an important consideration because conditional 
stability criteria typically link the maximum permissible time step to the local grid spacing rather 
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than to the underlying physical time and space "sampling" requirements of the solution. Stable 
implicit methods can benefit from choosing the time step to satisfy an accuracy requirement rather 
than a stability constraint, with the consequence that fewer time steps are needed for a fixed sim- 
ulation time. Accordingly, a one-dimensional model problem with nonuniform grid is developed 
to explore the influence of grid nonuniformity on solution accuracy and efficiency The explicit 
FDTD algorithm is used as a convenient reference algorithm for comparison. 

In the remainder of this report. Section 2 develops first and second order, characteristic based, 
finite difference LU/AF algorithms for Maxwell's equations in one dimension. The second or- 
der scheme includes the conduction current term. Section 3 gives a Fourier analysis of stability, 
damping and dispersion properties of the one dimensional algorithms. Section 4 discusses the 
characteristic based treatment of boundary conditions. Section 5 gives computed model problem 
results related to both accuracy and efficiency for propagation and scattering in free space and di- 
electric materials, using both uniform and nonuniform grids. Finally, some conclusions are drawn 
in Section 6. 


2 One-Dimensional LU/AF Algorithm 

This section details the theoretical background for development of both first and second-order 
accurate LU/AF algorithms. The second order algorithm is used in all comparisons with the 
FDTD method, and the first order algorithm is used at points adjacent to the outer boundaries of 
the computational grid. 


2.1 First-order algorithm 

Maxwell's equations for linear, homogeneous and lossless media in the one-dimensional case 
(taking d/dy = d/dz — 0) are 


dEy 

at 

dH z 


+ 


1 dHz 

dx 

1 dE, 


+ -- 


dt ' n dx 

These equations can be rewritten in flux conservative form as 


= 0 

( 1 ) 

= 0 

( 2 ) 




(3) 


where q = [E y , H Z ] T , f = [H z /e , E y /y] T and T denotes transpose. 
LU / AF algorithm, the flux conservative form of Maxwell's equations in 
ing form 


dq 

dt 


+4 


-0 


To develop the upwind 
(3) is recast in the follow- 

(4) 


where A is the Jacobian of / and is given by 


A = 


d l 

dq 


0 1/e 

1 //* o 


(5) 
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The matrix A has eigenvalues Ai^ = ±1/ y/JIe corresponding to right and left propagating waves 
with speeds ±c. The eigenvalue matrix is given by 


A = 


A l 0 

o a 2 


(6) 


The matrix A can be obtained from the eigenvalue matrix via a similarity transformation given by 

A = SAS- 1 (7) 

where the matrix S is composed of the eigenvectors of A and is given by 


5 = 


vWf -vWe 

5-1 _ 1 

i 

1 1 

’ 5 “2 

-s/tjji 1 


(8) 


The eigenvalue matrix A can be split into two parts, one each for the right and left propagating 
waves and is given by 


A = A+ +A' 


A + = 


' A| 

0 ‘ 

A — 

' 0 

0 

0 

0 

1 A — 

0 

A2 


Substituting (9) into (7) gives 

A = S {A + + A~) S~ l = A + + A- 
A ± = SA ± S ~ 1 

The flux vector / is then split into two parts given by 

/ = f + + r 

/ _± = A ± q 


( 9 ) 

( 10 ) 


( 11 ) 

(12) 


(13) 

(14) 


This flux-vector-splitting method is similar to that developed by Steger and Warming [28] for the 
Euler equations governing inviscid fluid flow. To construct the LU/AF algorithm, time and space 
are discretized to t n = nAt, = i Ax and the term Of /Ox is approximated at t = t n + (3 At by 


71-j-l 




( 


fi-fi-i I fJ±izJL 

Ax Ax 


(15) 


Here, [3 (0 < /3 < 1) is a parameter for evaluating the spatial approximation at t n + (3 At. This will 
give an explicit scheme for j3 = 0 and an implicit scheme for (3 > 0. The finite difference equation 
for (3) is given by 


At 


+ 0 


fi-fi-i , fr + 1 - fi 


p— \ n+l 


Ax 


+ 


Ax 


+ ^ + fl+1 — fl 


Ax 


Ax 


= 0 (16) 


Making the substitution Aqf = r/" +1 — qf , equation (16) can be rearranged for a uniform grid to 
give 


[//A* + f> (a^ : (A*-) + A+ (i--))] A,t = -til (17) 


4 



where / is the 2 by 2 identity matrix, and 


A h(‘) 
Ah(-) 


Ax 

.. (Oi+i - (•),- 


Ax 


( 18 ) 

(19) 


are the first-order backward and forward difference operators, respectively. The term is called 
the "residual" and is defined by 

Rl = ^J+' n + AtJ-' n ( 20 ) 

Note that equation (17) is an unfactored, implicit scheme of O [Ax, (ft — 1/2) At, At 2 ] accuracy for 
the numerical solution of Maxwell's equations. The solution of (17) involves solving a block tridi- 
agonal system of equations. The solution of tridiagonal systems can be avoided by an approximate 
factorization (AF) of the left side of (17) which results in the following LU/AF algorithm: 

[l/At + ftAft(A+.)}Aq* = -Rl (21) 

[I/At + 0A+ (A"-)] A qf = I/AtAq* (22) 

where A q* is an intermediate solution vector. This LU/AF algorithm is closely related to those 
developed previously for Computational Fluid Dynamics (e.g. [27]). It can be shown that these 
equations are a consistent and unconditionally stable approximation to (3), but the details of that 
analysis are omitted here for the sake of brevity. Note that the solution involves a two-step pro- 
cedure. In solving (21), a sweep through the grid in the forward direction results in a block lower 
bidiagonal matrix and a sweep through the grid in the backward direction results in a block upper 
bidiagonal matrix for (22). These equations are then solved by forward and backward substitu- 
tion, respectively. Thus, a costly inversion of a block tridiagonal matrix is avoided, and each step 
in the solution is effectively explicit. Substituting equation (22) into (21) and expanding gives, 

[;/Af + D (Af, (i+.) + A+ (i--)) + ft 2 At (a-, (A+-) A+ (yt-.))] Ag" = (23) 

Note that equation (17) is recovered except for the third term in brackets, which represents the fac- 
torization error. In this particular example, and A~ are constant matrices such that A + A~ = 0, 
and the factorization error is zero. However, the factorization error is nonzero in the following 
example with conduction currents, and more generally in two and three dimensional implemen- 
tations. This error can be reduced or eliminated by iterative error reduction. The present one- 
dimensional model problem results do not include factorization error except for the case including 
conduction currents. 


2.2 Second-order algorithm 

To develop the O (A t 2 , Ax 2 ) LU/AF algorithm, the flux conservative form of (4) is rewritten 
to include the electric and magnetic conduction currents as 


dq t dq 

it +A T x 


~ ~Pq 


(24) 
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where P is given by 


P = 


a/e 0 
0 a* /(i 


(25) 


and A is given by equation (5). The time derivative in (24) is approximated by a /3-weighted, 
0(At 2 ) difference equation given by 


Oq ~ (20 +^^-40^ + ( 2 ^-l) g f- 1 

Ot ~ 2 At [ ) 

The spatial derivative is again replaced by a /^-weighted average between time level n + 1 and n 
as in (15). This can be rewritten using the operator notation as 


(a 2 -£ + + AtJf) n+l + (1-0) (a nf? + a + f-y (27) 

The A 2t operators are now 0( Ax 2 ) difference operators to be defined shortly. The parameter 0 can 
be used to construct a series of explicit and implicit schemes. For example, if 0 = 0, this results 
in a leapfrog scheme; 0 = 0.5 results in a Crank-Nicolson scheme and 0=1 results in an Euler 
implicit scheme. Using (26) and (27), the finite-difference equation for (24) is 


(2/? + l)Agf J£0 1)A# + p ( A _j . f + A + f-y +l + { i-0) (A 2 - + A+f-y = 

-0P<% + '-(l-0)P% (28) 

Using the same flux vector splitting as in (14), this can be rearranged as 

[//(2Af) + 0aP + 0a (A 2i A+(-) + A+ A"(-))] Aq? = -R n 2l (29) 


The residual, R 2i , is now defined by 


R%i = c*\Pqi+ Ay+' n + a+/; 


-,n 2/3-1. „ ! 


2Af 


A 9 ? 


1 


(30) 


where a = 1/(20 -I- 1). The difference operators in equation (30) are replaced by 0(Ax 2 ) backward 
and forward upwind difference operators on three-point one-sided stencils defined by 


A 2 - (•) = (3(-)i - 4(-)»-i + (-)i-2) /(2Aar) (31) 

A+(.) = (-3(.)< + 4(-)i+i - (■)*+ 2 ) /(2Ax) (32) 


Substituting (14) into (30), the residual is given by 

K = n {rtf + A 2 -A + ?t + r 1 } (33) 

Equation (29) is an 0(At 2 , Ax 2 ), unfactored, upwind scheme for electromagnetics. The LU/AF 
scheme is defined by factoring the left side of (29) into two operators, each designed for a forward 
and backward grid sweep as in the first order implementation. The LU/AF scheme is then given 
by 

[l/(2At)+0aP + 0aA»A+(-)]Aq* = -Rg (34) 

[l/(2At) + 0aP + 0aAt l A-(-)]Aq? = [I / (2At) + 0aP] Aq* (35) 
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3 Fourier Analysis 

A Fourier analysis shows that both the first and second-order upwind LU/AF algorithms are 
unconditionally stable for /? > 1/2, and that they contain both numerical dissipation (or damping) 
and dispersion. The dissipation is present due to even order spatial derivatives in the truncation 
error which are a result of the upwind approximation. These schemes are consistent approxima- 
tions with 0((/3 - 1/2) At, A t 2 , Ax) truncation error for the first-order algorithm and 0(At 2 , Ax 2 ) 
for the second-order algorithm. A Fourier analysis is now given for the first-order LU/AF algo- 
rithm; the second-order algorithm analysis follows a similar development. 


Assuming an equally spaced mesh with periodic boundary conditions, a trial solution of the 
form 


qf = q be J '( n *-"> 

(36) 

is substituted into equations (21) and (22) to obtain: 


TAqf = Uq\ f 

(37) 

The matrices T and U are given by 


T = / + ^ (l - e>") I + A (e-jl ^ I'j A-] 

(38) 

v = -(2(i-^)^ + 2(.-i*_i)i-) 

(39) 

where c = 1/^/JIe, u is the Courant number defined by // — cAt/Ax and ()> = 
the substitution A q? = ^ n+1 — , equation (37) can be rearranged to give 

u6. Again making 

II 

i£r 

(40) 

where the amplification matrix G is defined by 


to 

7 

•Eh 

+ 

1^ 

II 

to 

(41) 


The stability of this scheme is governed by the spectral radius of the amplification matrix G, and 
its eigenvalues, G\ and Gz- Closed form solutions were obtained for G\ and Gz, but they are 
omitted here for the sake of brevity. A rigorous numerical stability analysis was also performed 
on the second-order LU/AF algorithm using the MAPLE 6™ software package. Checking the 
magnitude of the eigenvalues G \ and Gz for all values of /3, u and 0 showed that both the first 
and second-order schemes are unconditionally stable for ^ > 1/2. The stability analysis in this 
particular case is identical to the unfactored implicit scheme given in equation (17), because the 
factorization error is zero for P = 0. 



(42) 

(43) 
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Equation (43) is then solved for <f>*, the numerical wavenumber. The quantity 3m {<p*/<f>} gives 
the numerical dispersion, or phase velocity error. Extensive numerical experiments were per- 
formed to analyze the dispersive properties of these algorithms. We let 6 = 2n/N, where N is the 
grid resolution in cells/A. Figure 1 shows the numerical dispersion versus grid resolution with 
v = 3/4 for the first and second-order LU/AF algorithms as compared to the FDTD algorithm, 
which is 0(At 2 , Ax 2 ). The Courant number of v = 3/4 was chosen for these methods simply as 
a representative Courant number which could be expected for a problem involving a nonuniform 
grid. For a nonuniform mesh with a mesh stretch ratio of M s : 1 where M s = Ax max /Ax m i n/ 
the Courant number for FDTD varies pointwise in the grid over the range 1/M S < Vi < 1, where 
Vi = (cAt)/Axi. The numerical experiments showed that the first-order LU/AF algorithm has 

10 0 4 

o80: 


■H 

w 
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^2 0 


dP i . .. J 

10 20 30 40 50 60 70 80 

Grid Resolution { cells/lambda) 

Yee 

First Order 

Second Order 


Figure 1: Numerical dispersion of FDTD, first and second-order LU/AF algorithms versus grid 
resolution for /? = 1/2 and v = 3/4. 

lower dispersion errors than the second-order LU/AF method for v < 1. This is not particularly 
troublesome, because for v < 1, explicit schemes are generally more efficient and are preferred. 
However, for v > 1, explicit schemes cannot be used and the second-order LU/AF algorithm does 
have lower dispersion error than the first-order LU/AF method. For large values of v, the dis- 
persive error for the first and second-order LU/AF schemes begin to converge as shown in Figure 
2. The results also showed that the lowest numerical dispersion is obtained when /? = 1/2 as 
shown in Figure 3. Since the LU/AF method uses windward differencing, numerical dissipation 
(or damping) is present in the solution. The normalized dissipation is obtained as 'Re {</>*/</}. Fig- 
ures 4-6 show numerical dissipation for the first and second-order LU/AF schemes versus grid 
resolution, v and ft, respectively. Note that the second-order algorithm has much lower dissipation 
than the first order method. 
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First Order 

Second Order 

Figure 2: Numerical dispersion of first and second-order LU/AF algorithms versus Courant num- 
ber v at P = 1/2 and N = 10 cells/ A resolution. 



First Order 
Second Order 


Figure 3: Numerical dispersion of first and second-order LU / AF algorithms versus at v = 2 and 
N = 10 cells/A resolution. 
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Grid Resolution (cells /lambda) 


First Order 
Second Order 


Figure 4: Normalized dissipation of first and second-order LU/AF algorithms versus grid resolu- 
tion at j3 — 0.5 and v — 0.75. 



First Order 
Second Order 


Figure 5: Normalized dissipation of first and second-order LU/AF algorithms versus Courant 
number v at /? = 1/2 and N — 10 cells/ A resolution. 
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beta 

First Order 

Second Order 

Figure 6: Normalized dissipation of first and second-order LU/AF algorithms versus j3 at v = 2 
and TV = 10 cells/ A resolution. 

4 Boundary Conditions 

4.1 Outer Radiation Boundary Condition 

The characteristic based LU / AF method requires no extraneous boundary condition such as 
the Liao absorbing boundary condition [29] or the PML [30]. The FDTD method uses a spatial 
central difference operator, which for a wave propagating from left to right, eventually requires a 
grid point outside the domain. This requirement introduces an additional equation (i.e. boundary 
condition) to solve the system and introduces information into the solution that is not required 
by Maxwell's equations. Using an upwind characteristic based approach, the interior point algo- 
rithm calculates the left-going characteristic at the left boundary (i.e. at i = 0) and the right-going 
characteristic at the right boundary (i.e. at i = imax). Therefore, the only additional information 
required is information about waves that are entering the domain. Waves exiting the domain are 
naturally handled by the interior point algorithm. Therefore, the characteristic boundary condi- 
tions are implemented as follows: at grid point i = 0, equation (22) is used along with a specifica- 
tion of the incoming, right-going flux, /q\ At grid point i = imax, equation (21) is used along with 
a specification of the incoming, left-going flux, f~ max . Therefore, the only additional information 
introduced at the boundary is nothing more than what is required by the physical system. In mul- 
tidimensional problems, the local coordinates at the outer boundaries are rotated to align with the 
direction of wave propagation defined by E x H. The characteristic equations are developed along 
this direction and are appropriately applied at the boundaries. This procedure was discussed and 
outlined by Shang [11]. 
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4.2 Dielectric Surface Boundary Condition 

Since the LU/AF scheme follows the direction of information propagation (i.e. the character- 
istic), at a material interface, the slope of the characteristic curve (i.e. the speed of light in the 
material) changes. Therefore, for the LU/AF method to be widely applicable, a careful treatment 
of material interfaces is required. With a material boundary in place, the right-going and left-going 
characteristics see a change in characteristic speeds, and therefore, a material interface condition 
needs to be implemented to correctly model the physics. To implement this feature, consider the 
one-dimensional grid shown in Figure 7 where a dielectric boundary has been inserted at grid 
point if,. For the second order LU/AF method, at an arbitrary grid point i, the windward differ- 



Figure 7: One-dimensional FDTD grid showing a dielectric boundary at grid point if>. 

encing uses flux components / ; + 2 , /+_-[, /+, /“, f~ + 1 , ff +2 - This creates a difficulty at a dielectric 
interface since the algorithm assumes the material is homogeneous. We know at grid point if, that 
the physical boundary conditions require that 

Efan 1 

FI tan 1 

which requires that 

Ey\ 

H z i 

Imposing these boundary conditions at i i then makes the solution vector q constant between ma- 
terial 1 and 2. Therefore, the same solution vector q can be used for both materials at grid point 


= E t9n 2 (44) 

= H tan 2 (45) 

= E y 2 (46) 

= H z 2 (47) 
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ib. In general, in regions 1 and 2, we require flux components /+, /f , f£ and / 2 ~, which are re- 
stricted to the respective materials and are dependent on the material properties. However, we 
define additional flux components /f located at grid points ib, ib + 1 and ib + 2 in material 2 and 
flux components f.f at grid points ib and ib — 1 in material 1. By defining these components, we 
can then apply the second order windward difference formulas at grid points ib — 2, — 1, ib, 

ib + 1 and ib + 2 to simulate a homogeneous material region. For example, on the forward sweep, 
at grid point i b , we use the components /j + (i b - 2), /, + ( i b - 1), /+ {ib), /f {ib), /f (ib + 1)/ /f (*6 + 2) 
to update the solution vector q using the material properties of region 1. Similar analogies can be 
drawn at grid points *6 — 2, t& — 1, *& + 1 and ib + 2. These additional fluxes in each region in the 
vicinity of the dielectric boundary are computed from (14) based upon the material properties in 
each region, and they bridge the solution between regions 1 and 2. 

5 Model Problem Results 

The model problem results presented in this section demonstrate that the one dimensional im- 
plicit characteristic-based LU/AF algorithm can produce accurate results on a nonuniform grid 
for time steps significantly larger than the maximum permissible for a typical conditionally stable 
scheme. Although the operation count for the implicit scheme is higher, as the mesh variation is 
increased, fewer time steps are needed for an accurate solution over a fixed simulation time, and 
this savings more than offsets the additional operations (beyond a certain mesh variation). 

The second-order LU/AF algorithm was tested by implementing equations (34) and (35) for 
interior grid points away from absorbing boundaries and equations (21), (22) at grid points next 
to the absorbing boundaries. Characteristic-based boundary conditions were used to terminate 
the computational domain and the incoming flux (f^ max ) at the right boundary was set to zero. 
The code was initialized by writing a time snapshot of a propagating pulse in the grid. Several 
different types of problems were analyzed with the LU/AF algorithm and compared with the 
FDTD method in an attempt to assess the characteristics of the LU/AF algorithm as applied to a 
one-dimensional problem with a non-uniform grid. These are outlined in the following sections. 

5.1 Propagation 

To assess the dispersion and dissipation characteristics of the LU/AF algorithm, several propa- 
gation problems were analyzed using free space and dielectric materials and were compared with 
FDTD. A Gaussian pulse of the form 


E yt = e V ” > (48) 

was used as an excitation source for both propagation and scattering problems. For propagation, 
the code was changed to use periodic boundary conditions to enable simulation of propagation 
over large distances compared to the shortest wavelength contained in the frequency spectrum 
of the pulse. To illustrate the potential benefit of the LU/AF algorithm over conventional explicit 
schemes, the main parameters of interest are the time resolution (i.e. number of time steps /period) 
and the grid resolution (i.e. number of cells/wavelength) of the highest frequency of interest in 
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any given problem. These parameters will be designated as N t and N x , respectively. The other 
parameter of interest is the mesh stretch ratio, M s . The highest frequency of interest, f max , is 
calculated based upon N x and the largest cell size, and the time step for the implicit scheme is 
calculated based upon f max and the desired time resolution, N t . For the FDTD method, the time 
step was calculated based upon the Courant stability condition using the smallest grid size. 

The first propagation problem involves propagation in free space on a uniform mesh. The 
time and grid resolutions are N t = 30 and N x = 30, the cell size was 1 cm, the time step was 33.3 
ps and u — 1. The maximum frequency of interest was 1 GHz and the problem space was 2000 
cells. The pulse was allowed to propagate a distance equal to 100 wavelengths of the minimum 
wavelength, which corresponds to a distance of about 30 meters. Figure 8 shows the electric field 
versus distance obtained after 3000 time steps computed using the exact solution, FDTD and the 
LU/AF algorithm. Note the LU / AF scheme has slightly less accurate results due to the larger dis- 
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Figure 8: Electric field versus distance for free space propagation on a uniform mesh with (3 = 1/2 
and v = 3/4. 

persion error; but the general shape of the pulse is still acceptable. Again, this is not troublesome, 
because for v < 1, explicit schemes are more efficient and are preferred. Further numerical sim- 
ulations confirmed that the numerical dispersion was lower and the results were more accurate 
for the second-order LU/AF method for the case when v = 2 as shown in Figure 9. Additional 
simulations confirmed that the numerical dispersion error increases with Courant number v as 
with the ADI FDTD scheme [24]. 

To demonstrate the LU / AF scheme on nonuniform meshes, the same problem outlined above 
was simulated using mesh stretch ratios of M s =10:l and A/, -100:1 and the mesh was periodic ev- 
ery 10 cells. Figure 10 illustrates an expanded section of the nonuniform mesh for a mesh stretch 
ratio of M s =10:1. The smallest grid cells force the FDTD algorithm to take a time step 1 / 10th of 
the previous case, and it must run 10 times as many time steps to propagate the pulse the same 
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Figure 9: Electric field versus distance for free space propagation on a uniform mesh with 0 = 1/2 
and v = 2. 
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Figure 10: A section of a nonuniform mesh with a mesh stretch ratio of 10:1 and largest cell size of 

1 cm. 
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distance. However, the LU/AF algorithm can take the same time step as before, because the time 
resolution, N t — 30, has been preserved and the time step is set solely upon accuracy require- 
ments, not on the Courant stability condition. The maximum cell size remains constant, and since 
N x is also the same, the maximum frequency of interest remains constant at 1 GHz. The maximum 
grid coordinate changes from 20.0 m in the uniform mesh case to 8.17 m for the 10:1 nonuniform 
mesh. The Courant numbers in the grid for FDTD are in the range 0.1 < v < 1 and the Courant 
numbers for the LU/AF algorithm are in the range 1 < v < 10. Figure 11 shows the electric 
field versus distance obtained with the exact solution and after 30,000 time steps computed using 
FDTD and after 3,000 time steps using the LU/AF algorithm. Note the agreement is very good 
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Figure 11: Electric field versus distance for free space propagation on a nonuniform mesh with 
/3 = 1/2 and v = 1. 

for both methods in this case. In fact, Figure 12 shows the error in electric field versus distance 
for both methods and the errors are almost equal in magnitude. The LU/AF method still required 
about 15 seconds, but the FDTD method required about 22 seconds. Note the improved perfor- 
mance of the LU/AF algorithm away from the area of the main pulse. Since the mesh stretch ratio 
was 10 for this problem, the FDTD method is forced by the Courant stability condition to take 10 
times the number of time steps to simulate the same amount of physical time. It appears that the 
crossover point in efficiency between the LU/AF and FDTD methods is approximately a mesh 
stretch ratio of 10:1. This is in contrast to recent results for the ADI FDTD scheme [23] where the 
results suggested a crossover point of 30:1 mesh ratio. 

For a mesh stretch ratio of M s =100:1, the Courant number for FDTD was in the range 0.01 < 
v < 1 and for LU/AF was in the range 1 < u < 100. Execution time for FDTD was 230 seconds 
for 300,000 time steps. Execution time for the LU/AF method again required 15 seconds for 3,000 
times steps and the error in electric field for this simulation is shown in Figure 13. Note again 
that the LU/AF result appears "cleaner" than the FDTD result with less numerical noise. Based 
upon these results, we conclude that the LU/AF method is more efficient for either a.) nonuni- 
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Figure 12: Error in electric field versus distance for free space propagation on a nonuniform mesh 
with a mesh stretch ratio of 10:1 and (3 = 1/2, v — 1 (for FDTD). 
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Figure 13: Error in electric field versus distance for free space propagation on a nonuniform mesh 
with a mesh stretch ratio of 100:1 and /3 = 1/2, v = 1 (for FDTD). 
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form meshes with large mesh stretch ratios or b.) a uniform mesh with a low frequency incident 
excitation, where the pulse is highly resolved spatially over its entire frequency spectrum. Further 
computational experience on additional free space propagation problems on nonuniform grids re- 
veals that this sweeping method fails for large At (i.e. large Courant number). This can occur for 
two reasons, even though the Fourier stability analysis indicates unconditional stability. First, it is 
inapropriate to use a time step larger than one which propagates a wave over a distance compara- 
ble to the entire solution domain in a single time step, and this can cause instability. Secondly, the 
spatial sweeping process associated with each factor of the algorithm can itself become unstable 
if it loses diagonal dominance, as observed by Jameson and Turkel [25]. Although a first-order 
spatial upwind scheme guarantees diagonal dominance, there is a loss of diagonal dominance for 
higher-order upwind schemes as the time step becomes asymptotically large. The upwind scheme 
significantly increases the diagonal term compared with a centered scheme, however, and in ap- 
plying the LU/AF scheme, this problem can be avoided by exercising care in selecting the time 
step. 

5.2 Lossy Dielectric Materials 

To illustrate propagation in dielectric materials, the problem space was filled with a lossy di- 
electric material with parameters e = 4eo, M = Mo arid a = 0.02. The pulse was allowed to prop- 
agate on a uniform mesh for approximately 35 wavelengths at the minimum wavelength. Figure 
14 shows the results of FDTD versus LU / AF with excellent agreement between the two methods. 
A similar level of agreement is found with the same problem on a nonuniform grid. 



Figure 14: Electric field versus distance for propagation in a lossy dielectric on a uniform mesh 
with {3 = 1/2 and v = 1. 

To illustrate the use of the LU/AF scheme for scattering problems, reflection from a lossy 
dielectric half-space was considered as the benchmark problem. The dielectric interface scheme 
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of Section 4.2 was implemented and the results are again compared with the FDTD method. To 
solve this problem, the problem space is 2000 cells and it is filled with a lossy dielectric material 
from cell numbers 751-2000. The total electric field is recorded versus time at cell number 750. The 
incident field is obtained by running the code with free space only and recording the field versus 
time at the same location. The incident field is subtracted from the total field to give the scattered 
field. A point source located at the left boundary of the grid was used as the excitation source 
and the dielectric half-space material parameters are e — 4eo, n — no and a = 0.2 S/m. Figure 
15 shows the time-domain scattered field for this case, and the agreement is excellent. Figure 16 
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Figure 15: Scattered electric field versus time for scattering from a lossy dielectric half-space on a 
uniform mesh with j3 = 1/2 and v = 1. 

shows the reflection coefficient results and again the agreement is excellent. The same problem 
was run on the uniform mesh with a = 2, and on a nonuniform mesh with a mesh stretch ratio of 
10:1 and a = 0.2 S/m. Figures 17 and 18 show the respective results. 

Note that the LU/AF method produces more inaccurate results for moderate conductivity val- 
ues. This is due to the factorization error terms. As a is increased, the factorization errors become 
larger. In principle, this factorization error can be reduced throught iterative error reduction [31]. 
This has been demonstrated for CFD problems [27], but it is not explored in the present work. 
For perfect conductors, we simply set q — 0. Figure 18 shows that the LU/AF method is superior 
on nonuniform grids, even for problems involving lossy dielectrics. These results indicate that 
the physical dielectric boundary conditions have been properly implemented, thus making the 
LU/AF method applicable to a wider class of problems. 


6 Conclusions 

This report has introduced a LU/AF implicit finite-difference method for computational elec- 
tromagnetics. A second-order accurate, LU/AF, characteristic based algorithm for electromag- 
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Figure 16: Reflection coefficient magnitude versus frequency for scattering from a lossy dielectric 
half-space on a uniform mesh with /? = 1/2 and v = 1. 
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Figure 17: Reflection coefficient magnitude versus frequency for scattering from a lossy dielectric 
half-space on a uniform mesh with (i = 1/2, v = 1 and cr = 2 S/m. 
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Figure 18: Reflection coefficient magnitude versus frequency for scattering from a lossy dielectric 
half-space on a uniform mesh with 0 = 1/2, v = 1 and a — 0.2 S/m. 


netics has been implemented and tested on one-dimensional model problems for uniform and 
nonuniform grids. The one-dimensional model problem results for the characteristic based im- 
plicit scheme demonstrate that: 

1. Accurate solutions for wave propagation can be obtained using Courant number signifi- 
cantly greater than one on nonuniform grids. 

2. Stability of the implicit scheme allows a given amount of time to be covered in fewer time 
steps than a conditionally stable scheme, and this can significantly improve efficiency when 
nonuniform grids are needed. This is also extremely beneficial for low frequency excitation 
sources where geometric details may force a highly oversampled mesh (uniform or nonuni- 
form) in certain regions. 

3. The method works well with lossy dielectric materials and perfect conductors. 

4. The factorization error increases as a increases, but it can be reduced in principle through 
iterative error reduction. 

5. The lowest dispersion and damping errors occur when 0 = 1/2. 

The present results demonstrate potential advantages in a one dimensional context, and this 
approach appears promising for development of stable, accurate and efficient implicit LU/AF 
schemes for complex two and three dimensional applications. Extensions to two and three-dimensional 
applications have been outlined [31], but details of this work will the subject of future reports and 
articles. 
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